How to code trigonometry question to solve Milne-Simpson method?

[arin mn]

[arin mn]

clear;clc;

syms('t','y');

s=('k1*y.^2+K2');

f=inline(s,'t','y');

N=10;

S=10;

t=zeros(1,N+1);

y=zeros(1,N+1);

r=zeros(1,S+1);

OUP=1;

a=0;b=1;

t(1)=0;%for t=0.0

r(1)=0;%for r=0.0

R=0.1; %step size of r

h=(b-a)/N;

fprintf(OUP,'step size(h)=');

disp(h);

fprintf(OUP,'\nThe exact solution for min of function is = (0.25+0.75r)e^t\n');

fprintf(OUP,'The exact solution for max of function is =(1.125-0.125r)e^t\n');

for j=1:S+1 %for looping of r : 0<r<1

Y0=(0.5+0.5*r(j)); %initial condition of y min

Z0=(0.0012-0.0012*r(j)); % initial condition of ymax

%max

y(1)=Y0;

z(1)=Z0;

%%EXACT SOLUTION

yExact=(d*r(j))*tan(e(t));

zExact=(f*r(j))*tan(g(t));

fprintf(OUP,'\n\nAt R= %2.1f \n',r(j));

fprintf(OUP,'\nEXACTSOLUTION\n');

fprintf(OUP,'______________________________________________________\n');

fprintf(OUP,' step t y(min) y(max)\n');

fprintf(OUP,'______________________________________________________\n');

fprintf(OUP,'%3d %11.2f %11.15f %11.15f\n',1,t(1),yExact(1),zExact(1));

for i=1:N

t(i+1)=t(i)+h;

yE(i+1)= (d*r(j))*tan(e(t+1));

zE(i+1)=(f*r(j))*tan(g(t+1));

fprintf(OUP,'%3d %11.2f %11.15f %11.15f\n',i+1,t(i+1),yE(i+1),zE(i+1));

end

%APPROXIMATE SOLUTION

fprintf(OUP,'\nMILNE SIMPSON METHOD\n');

fprintf(OUP,'______________________________________________________\n');

fprintf(OUP,' step t y(min) y(max)\n');

fprintf(OUP,'______________________________________________________\n');

fprintf(OUP,'%3d %11.2f %11.15f%11.15f\n',1,t(1),y(1),z(1));

%Runge-Kutta

for i=1:3

t(i+1)=t(i)+h;

%formula Runge-Kutta

%for ymin

k1=h*f(t(i),y(i));

k2=h*f(t(i)+h/2, y(i)+(k1/2));

k3=h*f(t(i)+h/2, y(i)+(k2/2));

k4=h*f(t(i)+h,y(i)+k3);

y(i+1)=y(i)+(1/6*(k1+2*k2+2*k3+k4));

%for ymax

l1=h*f(t(i),z(i));

l2=h*f(t(i)+h/2, z(i)+(l1/2));

l3=h*f(t(i)+h/2, z(i)+(l2/2));

l4=h*f(t(i)+h,z(i)+l3);

z(i+1)=z(i)+(1/6*(l1+2*l2+2*l3+l4));

%fprintf(OUP,'%3d %11.2f %11.15f %11.15f\n',i,t(i),y(i),z(i));

fprintf(OUP,'%3d %11.2f %11.15f %11.15f\n',i+1,t(i+1),y(i+1),z(i+1));

end

%Milne-Simpson Method

for i=4:N

t(i+1)=t(i)+h;

%formula predictor: Milne

yP(i+1)=y(i-3)+(4/3)*h*(2*f(t(i-2),y(i-2))-f(t(i-1),y(i-1))+2*f(t(i),y(i)));

zP(i+1)=z(i-3)+(4/3)*h*(2*f(t(i-2),z(i-2))-f(t(i-1),z(i-1))+2*f(t(i),z(i)));

%formula corrector: Milne

y(i+1)=y(i-1)+(h/3)*(f(t(i-1),y(i-1))+4*f(t(i),y(i))+f(t(i+1),yP(i+1)));

z(i+1)=z(i-1)+(h/3)*(f(t(i-1),z(i-1))+4*f(t(i),z(i))+f(t(i+1),zP(i+1)));

fprintf(OUP,'%3d %11.2f %11.15f %11.15f\n',i+1,t(i+1),y(i+1),z(i+1));

end

%ABSOLUTE ERROR

fprintf(OUP,'\nABSOLUTEERROR\n');

fprintf(OUP,'__________________________________________\n');

fprintf(OUP,' step t y(min) y(max)\n');

fprintf(OUP,'_________________________________________\n');

yError=abs(yExact(1)-y(1));

zError=abs(zExact(1)-z(1));

fprintf(OUP,'%3d %11.2f %11.5e %11.5e\n',1,t(1),yError,zError);

for i=1:N

t(i+1)=t(i)+h;

yError=abs(yE(i+1)-y(i+1));

zError=abs(zE(i+1)-z(i+1));

fprintf(OUP,'%3d %11.2f %11.5e %11.5e\n',i+1,t(i+1),yError,zError);

end

r(j+1)=r(j)+R;

end