function y=f73FE(dt)

%solve dudt=A cos(2pi ft)^2 - ku^2

%parameters for problem

A=0.5; f=1; k=5; u0=0;

tend=10;

%problem set up

t=0:dt:tend;

Lt=length(t);

u=zeros(Lt,1);

u(1)=u0; %initial condition

for n=1:Lt-1

u(n+1)=u(n)+dt*(A*cos(2*pi*f*(t(n)))^2-k*(u(n))^2);

end

y=[t' u];

function y=f73RK(dt)

%solve dudt=A cos(2pi ft)^2 - ku^2

%parameters for problem

A=0.5; f=1; k=5; u0=0;

tend=10;

%problem set up

t=0:dt:tend;

Lt=length(t);

u=zeros(Lt,1);

u(1)=u0; %initial condition

for n=1:Lt-1

k1=dt*(A*cos(2*pi*f*(t(n)))^2-k*(u(n))^2);

k2=dt*(A*cos(2*pi*f*(t(n)+dt/2))^2-k*(u(n)+k1/2)^2);

k3=dt*(A*cos(2*pi*f*(t(n)+dt/2))^2-k*(u(n)+k2/2)^2);

k4=dt*(A*cos(2*pi*f*(t(n+1)))^2-k*(u(n)+k3)^2);

u(n+1)=u(n)+(k1+2*k2+2*k3+k4)/6;

end

y=[t' u];

function y=fdudt(dt,tend,u0)

close all;

%solve dudt=-u

t=0:dt:tend+dt;

Lt=length(t);

u=zeros(Lt,1);

u(1)=u0; %initial condition

tic;

for n=1:Lt-1

u(n+1)=u(n)+dt*(-u(n));

end

toc;

y=[t' u];

tex=0:0.01:tend;

uex=exp(-tex);

plot(t,u,tex,uex); legend('Forward Euler','exact')

figure; plot(t,exp(-t)-u');

function y=fdudtRK(dt,tend,u0)

%solve dudt=-u

close all;

t=0:dt:tend+dt;

Lt=length(t); u=zeros(Lt,1); u(1)=u0; %initial condition

for n=1:Lt-1

k1=dt*(-u(n));

k2=dt*(-(u(n)+k1/2));

k3=dt*(-(u(n)+k2/2));

k4=dt*(-(u(n)+k3));

u(n+1)=u(n)+(k1+2*k2+2*k3+k4)/6;

end

y=[t' u];

tex=0:0.01:tend;

uex=exp(-tex);

plot(t,u,'bd',tex,uex,'g'); legend('RK','exact')

figure; plot(t,exp(-t)-u','r+');

function dXdt=freac1(t,X)

k1=100; k2=10; k3=0.1;

dXdt=[-k1*X(1)+k2*X(2);k1*X(1)-k2*X(2)-k3*X(2); k3*X(2)];

function dXdt=freac2(t,X,T)

R=8.3145; %universal gas constant

Ea1=50000; Ea2=Ea1; Ea3=25000; %activiation energy for A, B, C

A1=8.2e10; A2=8.2e9; A3=2409; %pre exponential factors

k1=A1*exp(-Ea1/(R*(T+273)));

k2=A2*exp(-Ea2/(R*(T+273)));

k3=A3*exp(-Ea3/(R*(T+273)));

dXdt=[-k1*X(1)+k2*X(2);k1*X(1)-k2*X(2)-k3*X(2); k3*X(2)];

function dXdt=freac3(t,X)

R=8.3145; %universal gas constant

Ea1=50000; Ea2=Ea1; Ea3=25000; %activiation energy for A, B, C

A1=8.2e10; A2=8.2e9; A3=2409;

k1=A1*exp(-Ea1/(R*(20+8*t+273)));

k2=A2*exp(-Ea2/(R*(20+8*t+273)));

k3=A3*exp(-Ea3/(R*(20+8*t+273)));

dXdt=[-k1*X(1)+k2*X(2);k1*X(1)-k2*X(2)-k3*X(2); k3*X(2)];

function loop(n)

close all

u=0;

for i=1:n-1

u(i+1)=u(i)+1;

end

plot(u)

function loopm(n)

close all

u=zeros(n,1);

for i=1:n-1

u(i+1)=u(i)+1;

end

plot(u)

function dy2=noteasy(t,u)

A=0.5; f=1; k=0.5;

dy2=[A*cos(2*pi*f*t)^2 - k*(u^2)];

function d2x=pend(t,x)

g=9.81; L=0.5;

d2x=[x(2);-(g/L)*x(1)];

function d2x=pendnl(t,x)

global g L

d2x=[x(2);-g*L*sin(x(1))];

%m-file to solve reaction scheme

clear; close all

%Initial conditions

t0=0; tfinal=10; tchange=5;

A0=10; B0=100; C0=0;

X0=[A0;B0;C0];

temp1=20;

temp2=50;

%we need matlab to ignore the options set

%we do this by passing a blank matrix into matlab

opts=[ ];

%solve

[t1,X1]=ode15s(@freac2,[t0 tchange],X0,opts,temp1);

newX0=X1(end,:);

[t2,X2]=ode15s(@freac2,[tchange tfinal],newX0,opts,temp2);

t=[t1;t2];

X=[X1;X2];

plot(t,X); legend('A','B','C'); xlabel('Time'); ylabel('conc')

%m-file to call f73FE

clear; close all

dt=input('enter a value of dt: ')

tic; y1=f73FE(dt); toc;

tic; y2=f73FE(dt/2); toc;

subplot(2,1,1); plot(y1(:,1),y1(:,2),'r',y2(:,1),y2(:,2),'k'); legend('dt','dt/2');

xlabel('time'); ylabel('Solution');

ycomp=interp1(y2(:,1),y2(:,2),y1(:,1));

subplot(2,1,2); plot(y1(:,1),ycomp-y1(:,2)); legend('difference')

solf73RK

%m-file to call f73RK

clear; close all

dt=input('enter a value of dt: ')

tic; y1=f73RK(dt); toc;

tic; y2=f73RK(dt/2); toc;

subplot(2,1,1); plot(y1(:,1),y1(:,2),'r',y2(:,1),y2(:,2),'k'); legend('dt','dt/2');

xlabel('time'); ylabel('Solution');

ycomp=interp1(y2(:,1),y2(:,2),y1(:,1));

subplot(2,1,2); plot(y1(:,1),ycomp-y1(:,2)); legend('difference')

solreac2

%m-file to solve reaction scheme

clear; close all

%Initial conditions

t0=0; tfinal=10;

A0=10; B0=100; C0=0;

X0=[A0;B0;C0];

temp=input('enter reaction temperature in degree C: ');

%we need matlab to ignore the options set

%we do this by passing a blank matrix into matlab

opts=[ ];

%solve

[t,X]=ode15s(@freac2,[t0 tfinal],X0,opts,temp);

plot(t,X); legend('A','B','C'); xlabel('Time'); ylabel('conc')

solreac3

%m-file to solve reaction scheme

clear; close all

%Initial conditions

t0=0; tfinal=10;

A0=10; B0=100; C0=0;

X0=[A0;B0;C0];

%solve

tic;[t,X]=ode45(@freac3,[t0 tfinal],X0);toc

plot(t,X); legend('A','B','C'); xlabel('Time'); ylabel('conc')

%m-file to solve equation noteasy

clear; close all

tfinal=10; tstart=0;

ystart=0;

tic;[T,Y]=ode45(@noteasy,[tstart tfinal],ystart);toc

figure; subplot(2,1,1); plot(T,Y); title('Solution');

subplot(2,1,2); plot(T(2:end),diff(T)); title('time step used')

%m-file to solve simple pendulum

clear; close all

g=9.81; L=0.5;

%define ICs

t0=0; tfinal=10;

x0=0.1; dx0=0;

ic(1)=x0; ic(2)=dx0;

%solve

[t,y]=ode45(@pend,[t0 tfinal],ic');

x=y(:,1); %the solution to x is in the first column of y

dx=y(:,2); %the solution to dx is in the second column

%exact solution

tex=t0:0.01:tfinal;

xex=x0*cos(sqrt(g/L)*tex);

%plot

subplot(2,1,1); plot(t,x,'b',tex,xex,'r')

legend('numeric','exact');

xlabel('time'); ylabel('Angle (radians)')

xexn=x0*cos(sqrt(g*L)*t);

subplot(2,1,2); plot(t,xexn-x);

legend('difference'); xlabel('time')

%m-file to solve simple pendulum

clear; close all

global g L

g=9.81; L=0.5;

%define ICs

t0=0; tfinal=10;

x0=input('please enter the initial value of x (eg 0.1): ');

if isempty(x0)==1

disp('using defualt value 0.1');

x0=0.1;

end

dx0=0;

ic(1)=x0; ic(2)=dx0;

%solve

[t,y]=ode45(@pendnl,[t0 tfinal],ic');

x=y(:,1); %the solution to x is in the first column of y

dx=y(:,2); %the solution to dx is in the second column

%exact solution

tex=t0:0.01:tfinal;

xex=x0*cos(sqrt(g/L)*tex);

%plot

subplot(2,1,1); plot(t,x,'b',tex,xex,'r')

legend('numeric','simplified model exact');

xlabel('time'); ylabel('Angle (radians)')

xexn=x0*cos(sqrt(g/L)*t);

subplot(2,1,2); plot(t,xexn-x);

legend('difference'); xlabel('time')