function hw4solutions
clear all
%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')
disp('press any key to continue')
pause

%lorenz
global rho
rho=28;
tr=[0 30]; %seconds
init=[27 9 12]; %300 rabbits; 150 foxes
[t,y]=ode45(@rklorenz, tr, init);
plot(t,y(:,1),'r',t,y(:,2),t,y(:,3),'g')
xlabel('time')
legend('x','y','z')
figure
plot3(y(:,1),y(:,2),y(:,3))
disp('press any key to continue')
pause

%lorenz
global rho
rho=5;
tr=[0 30]; %seconds
init=[4 3 4]; %300 rabbits; 150 foxes
[t,y]=ode45(@rklorenz, tr, init);
plot(t,y(:,1),'r',t,y(:,2),t,y(:,3),'g')
xlabel('time')
legend('x','y','z')
figure
plot3(y(:,1),y(:,2),y(:,3))
pause

% problem 4
xlow=1;
xhigh=2;
solinit = bvpinit(linspace(xlow,xhigh,20),[1 0]);
sol = bvp4c(@bvpode,@bvpbc,solinit);
xint = linspace(xlow,xhigh);
Sxint = deval(sol,xint);
yanalyt=sin(xint).^2;
plot(xint,Sxint(1,:),xint,yanalyt,'o')
xlabel('x')
ylabel('y')

function p=rkfox(t,w)
alpha=0.01;
R=400;
p=zeros(2,1);
p(1)=2*w(1)*(1-w(1)/R)-alpha*w(1)*w(2);
p(2)= -w(2)+alpha*w(1)*w(2);
%
function p=rklorenz(t,w)
global rho
beta=8/3;
sigma=10;
x=w(1); y=w(2); z=w(3);
p=zeros(3,1);
p(1)=-beta*x+y*z;
p(2)= -sigma*(y-z);
p(3)= -y*x+rho*y-z;
% -----------------------------------------------
function dydx = bvpode(x,y)
dydx = [  y(2)
         2-4*y(1).^2/sin(x).^2 ];
% -----------------------------------------------
function res = bvpbc(ya,yb)
res = [  ya(2)-2*sin(1)*cos(1)
         yb(1)-sin(2)^2 ];