First of all, the `NDSolve` solution can be further improved:

    k = 5.;
    tf = 3.; 
    c[θ1_?NumericQ, θ2_?NumericQ] := 
      NIntegrate[(Sin[θ1 + t]^2. Cos[θ2 + t])/E^t, {t, 0., tf}, 
       Method -> {Automatic, SymbolicProcessing -> 0}];
    
    sol = NDSolve[{Derivative[1][θ1][t] == θ2[t]/E^(Sin[t]/10.) + c[θ1[t], θ2[t]], 
         Derivative[1][θ2][t] == (-k Sin[θ1[t]])/E^(Sin[t]/10.), 
         θ1[0.] == 1., θ2[0.] == 1.}, {θ1, θ2}, {t, 0., tf}]; // AbsoluteTiming
    (* {0.0730725, Null} *)

    ref = sol[[1, All, -1]] // ListLinePlot

Then, as to the implementation of collocation method, `NMinimize` is undoubtedly a bad choice. Just use `FindRoot`. I've also turned to Chebyshev–Gauss–Lobatto grid to improve the accuracy of the solution:

    CGLGrid[{xl_, xr_}, n_] := 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)]);
      
    npoints = 11;
    tpoints = CGLGrid[{0., tf}, npoints];
    myd = dMatrixLagrange[tpoints];

    SetAttributes[c, Listable];
    θ1init = {1.}; θ2init = {1.};
    tst = FindRoot[{arg1, arg2} |-> 
                     With[{t = tpoints, θ1 = θ1init~Join~arg1, θ2 = θ2init~Join~arg2}, 
                          Rest /@ {myd . θ1 - (θ2/E^(Sin[t]/10.) + c[θ1, θ2]), 
                                   myd . θ2 - (-k Sin[θ1])/E^(Sin[t]/10.)}], 
                   {#, #} &@{ConstantArray[1., npoints - 1]}]; // AbsoluteTiming
    (* {0.864978, Null} *)

    ListPlot[{tpoints, #}\[Transpose] & /@ Join[{θ1init, θ2init}, tst, 2], 
      PlotRange -> All]~Show~ref

![Mathematica graphics](https://i.sstatic.net/VDbAL.png)

Remark
-

1. I've made use of the [hidden syntax of `FindRoot`](https://mathematica.stackexchange.com/a/163273/1871), a more basic approach is:
   ```
     help[arg1_List, arg2_] := 
      With[{θ1 = {1.}~Join~arg1, θ2 = {1.}~Join~arg2}, 
       Rest /@ With[{t = tpoints}, 
         {myd . θ1 - (θ2/E^(Sin[t]/10.) + c[θ1, θ2]), 
          myd . θ2 - (-k Sin[θ1])/E^(Sin[t]/10.)}]]
    
     tst = {arg1, arg2} /. 
      FindRoot[help[arg1, arg2], {#, ConstantArray[1., npoints - 1]} & /@ {arg1, arg2}]
   ```

2. `` NDSolve`FiniteDifferenceDerivative `` is a better choice to build `myd` e.g.:

       myd = NDSolve`FiniteDifferenceDerivative[1, tpoints, 
         DifferenceOrder -> 2]["DifferentiationMatrix"]

       NDSolve`FiniteDifferenceDerivative[1, tpoints, 
            DifferenceOrder -> "Pseudospectral"]["DifferentiationMatrix"] - 
          dMatrixLagrange[tpoints] // Abs // Max
       (* 2.66454*10^-14 *)

3. `|->` is introduced in _v12.2_, if you're not yet in _v12.2_, use `\[Function]` instead.