function odepractice
clear all
%insect problem
tspan=[0 50];
wnot=10000;
[t,y]=ode45(@rkinsect,tspan,wnot);
plot(t,y(:,1)) 
xlabel('time')
ylabel('N')
disp('press any key to continue')
pause

%tank draining problem
tspan=[0 2480];
wnot=9;
[t,y]=ode45(@rktank,tspan,wnot);
plot(t,y)
xlabel('time (s)')
ylabel('h (ft)')
disp('press any key to continue')
pause

%rocket problem
tspan=[0 40/0.8];
wnot=0;
[t,y]=ode45(@rkrocket,tspan,wnot);
plot(t,y)
xlabel('time (s)')
ylabel('V (m/s)')
disp('press any key to continue')
pause

%pendulum problem
tr=[0 5]; %seconds
init=[pi/8 0]; %start at pi/8
[t,y]=ode45(@rkpend, tr, init);
plot(t,y(:,1))
ylabel('theta')
xlabel('time(s)')
hold on
init=[pi/2 0]; %start at pi/2
[t,y]=ode45(@rkpend, tr, init);
plot(t,y(:,1))
init=[pi-0.1 0]; %start at pi-0.1
[t,y]=ode45(@rkpend, tr, init);
plot(t,y(:,1))
hold off
disp('press any key to continue')
pause

%fox/rabbit problem
tr=[0 15]; %seconds
init=[300 150]; %300 rabbits; 150 foxes
[t,y]=ode45(@rkfox, tr, init);
plot(t,y(:,1),'r',t,y(:,2))
ylabel('number of animals')
xlabel('time')
legend('rabbits','foxes')

function p=rkinsect(t,w)
R=0.55;
C=1e4;
nc=1e4;
r=10000;
p= R*w(1)*(1-w(1)/C)-r*w(1)^2/(nc^2+w(1)^2);

function p=rktank(t,h)
p= -0.0334*sqrt(h)./(10*h-h.^2);

function p=rkrocket(t,w)
T=48000;
g=9.81;
r=0.8;
b=40;
mo=2200;
m=mo*(1-r*t/b);
p= T/m-g;

function p=rkpend(t,w)
r=1;
g=9.81;
p=zeros(2,1);
p(1)=w(2);
p(2)= -g/r*sin(w(1));

function p=rkfox(t,w)
alpha=0.01;
p=zeros(2,1);
p(1)=2*w(1)-alpha*w(1)*w(2);
p(2)= -w(2)+alpha*w(1)*w(2);