In[]:=
Clear[U]
In[]:=
n=2;
In[]:=
vecU=Table[U[i],{i,0,n}];
In[]:=
ODEeqn
0
=
1
2
(1+4
2
t
)
(-243.5vecU[[3]][t]+363.19vecU[[2]][t])+
′′
vecU[[1]]
[t];​​
ODEeqn
1
=
1
2
(1+4
2
t
)
(-43.25vecU[[1]][t]-3.219vecU[[3]][t])+
′′
vecU[[2]]
[t];​​
ODEeqn
2
=
1
2
(1+4
2
t
)
(16.5vecU[[2]][t]+18.45vecU[[1]][t])+
′′
vecU[[3]]
[t];
In[]:=
ODETable=Table[
ODEeqn
i
0,{i,0,n}];(*Print[ODETable];*)
In[]:=
step=1.0;
In[]:=
t0=0;​​tf=15;​​IClist1=Table[U[i][0]0.01,{i,0,n}];IClist2=Table[U[i]'[0]0,{i,0,n}];​​IClist=Join[IClist1,IClist2];​​ODEandICTable=Join[ODETable,IClist];Vars=Table[U[i],{i,0,n}];​​currentSolU=AbsoluteTiming[NDSolve[ODEandICTable,(*WhenEvent[Mod[t,step]0,Print["t =",t]],*)Vars,{t,t0,tf},AccuracyGoalMachinePrecision/2,PrecisionGoal1+MachinePrecision/2,InterpolationOrderAll,MaxSteps1000000,CompiledFalse]];​​Print["Final Time = ",tf," ","Comp. time = ",currentSolU[[1]]]
Final Time = 15 Comp. time = 0.003651
In[]:=
For[i=0,i<n+1,i++,U[i][t_]=U[i][t]/.currentSolU[[2]]];
In[]:=
Plot[{U[0][t],U[1][t],U[2][t]},{t,0,tf},PlotRangeAll,AxesLabel{"t","
u
j
​[t]"},PlotStyle->{Black,Blue,Green}]