% Declare the relevant parameters
h=0.01;
w0=[0.0 0.0 0.0 1.0];
%w0=[0.0 0.000000001 0.0 100000000.0];
beta=2.0;
xmin=0.0;
xmax=1.0;

nmax=(xmax-xmin)/h;

w_euler_last = w0;
output = [0.0 w0];

% Beginning of Loop

for i=2:(nmax+1)

% Euler Method

    x = xmin + (i-1)*h;
    w_euler = w_euler_last + h*deriv(x,w_euler_last);
    w_euler_last =  w_euler;


output = [output; x w_euler];
end 

output2 = [0.0 w0 0.0 0.0];
for i=2:(nmax+1);

y = output(i,2) + (beta - output(nmax+1,2))/output(nmax+1,4) * output(i,4);
x = xmin + (i-1)*h;
yexact = exact(x);

output2 = [output2; output(i,1)  output(i,2)  output(i,3)  output(i,4)  output(i,5) y yexact];

end

% Plot Results for w, y

plot(output2(:,1),output2(:,6),'r',output2(:,1),output2(:,7),'b');

% Label the plot 

title('Linear Shooting: y''=4*(y-x)  with y(0)=0, y(1)=2');
xlabel('Position');
ylabel('y');
legend('Linear Shooting','Exact')
