Return to Revisions

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

Remark

1. I've made use of the hidden syntax of FindRoot, 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.