%code 2.2 for lecture 2; need u_0(x), u_L(t), u_R(t) dt=input(' dt= '); T=input( 'T= '); epp=input( 'epp= '); x1=1; x2=-1; for N=3:10 j=[1:1:N-1]; x=[cos(j*pi/N)]'; U=epp*(x+tan(0.5*x)); D1=DM1(N); D2=D1*D1; A=D2(2:N,2:N); B=D1(2:N,2:N); time=0.0; while time < T U0=U; C=diag(U); guo=1+epp*time; UR=(epp/guo)*(x1+tan(0.5*x1/guo)); UL=(epp/guo)*(x2+tan(0.5*x2/guo)); b=RK_RHSb(N,UL,UR,epp,D1,D2,U); K1=epp*A*U-C*B*U+b'; U=U0+0.5*dt*K1; C=diag(U); guo=1+epp*(time+0.5*dt); UR=(epp/guo)*(x1+tan(0.5*x1/guo)); UL=(epp/guo)*(x2+tan(0.5*x2/guo)); b=RK_RHSb(N,UL,UR,epp,D1,D2,U); K2=epp*A*U-C*B*U+b'; U=U0+0.5*dt*K2; C=diag(U); b=RK_RHSb(N,UL,UR,epp,D1,D2,U); K3=epp*A*U-C*B*U+b'; U=U0+dt*K3; C=diag(U); guo=1+epp*(time+dt); UR=(epp/guo)*(x1+tan(0.5*x1/guo)); UL=(epp/guo)*(x2+tan(0.5*x2/guo)); b=RK_RHSb(N,UL,UR,epp,D1,D2,U); K4=epp*A*U-C*B*U+b'; U=U0+dt*(K1+2*K2+2*K3+K4)/6; time=time+dt; end %exact solution guo=1+epp*T; exact=(epp/guo)*(x+tan(0.5*x/guo)); error=abs(exact-U); xx(N-2)=N; err(N-2)=max(error); ss=max(error) end fprintf(1, '%9.0f %10.2e \n', [xx; err]);