请问我想做下图这种图，但是代码报错不知道什么原因，求大佬指点一下。感谢
代码问题如下：

不好意思没看到发帖注意事项，一下为代码
GetR[k_]:=Module[{temp,j,n,sint},
sint = (Exp[omega*t*I]-Exp[-omega*t*I])/2/I;
temp[1] = Sum[a[n]*aa[k-1-n],{n,0,k-1}];
temp[2] = at[k-1]-Q1*(1+delta*sint)*a[k-1]+Q2*temp[1];
R[k] = Expand[temp[2]];
];
a1 = Sqrt[Q1/Q2];
a[0] = a1 + (a0-a1)*Exp[-gamma*t];
A[0] = a[0];
At[0] = D[U[0],t];
chi[k_]:=If[k<=1,0,1];
GetAll[k_]:=Module[{temp},
at[k] = D[a[k],t]//Expand;
aa[k] = Sum[a[n]*a[k-n],{n,0,k}]//Expand;
];
L[f_]:= D[f,t] + gamma*f;
Linv[f_]:=Block[{temp,g,G,s},
temp = DSolve[L[g[t]] == f,g[t],t]/.C[_]->0;
G = g[t]/.temp;
s = G[[1]]//Expand;
Expand[s]
];
Linv[p_Plus] := Map[Linv,p];
Linv[c_*f_] := c*Linv[f] /; FreeQ[c,t];
GetErr[k_,tmax_]:=Module[{temp,sint},
sint = (Exp[omega*t*I]-Exp[-omega*t*I])/2/I;
error[k] = At[k] - Q1*(1+delta*sint)*A[k] + Q2*A[k]^3//Expand;
err[k] = Integrate[error[k]^2,{t,0,tmax}]/tmax;
];
ham[m0_,m1_]:=Module[{temp,k,n,C1},
For[k=Max[1,m0],k<=m1,k=k+1,
Print[" k = ",k];
GetAll[k-1];
GetR[k];
RHS[k] = hbar*R[k]//Expand;
Special = Linv[RHS[k]];
temp = Special + chi[k]*a[k-1];
C1 = -temp/.t->0;
a[k] = temp + C1*Exp[-gamma*t]//Expand;
A[k] = A[k-1] + a[k]//Expand;
At[k] = D[A[k],t];
If[PRN == 1,
GetErr[k,tmax];
Print[" error = ", err[k]//N];
];
];
Print["Successful !"];
];
exact = 1/Sqrt[Q2/Q1-(Q2/Q1-1/a0^2)*Exp[-2*Q1*t]];
(* Input parameters *)
a0 = 1;
Q1 = 1;
Q2 = 4;
P. G. SIDDHESHWAR
Copyright © 2010 SciRes. AM
554
delta = 1/5;
omega = 10;
hbar = -1/2;
gamma = 2*Q1;
PRN = 0;
tmax = 10;
(* Print out the input parameters *)
Print[" a0 = ",a0];
Print[" Q1 = ",Q1];
Print[" Q2 = ",Q2];
Print[" delta = ",delta];
Print[" omega = ",omega];
Print[" hbar = ",hbar];
(* get fifth-order approximation *)
ham[1,5];

注意那个乘号×

EXP无定义，是Exp的笔误吧。
应该还有更严重的问题，不过一眼暂时看不出。