%Solve BVP by Matlab function and plot solution y and y'

function [] = bvp1
    solinit = bvpinit(linspace(0,1,5),[2 2]); %generate initial solution 
    sol = bvp4c(@bvp_def,@bvp_bc_def,solinit);

    figure(1);
    x = linspace(0,1); % artificial *time*, to show what is going on in one dimension.
    y1 = deval(sol,x); % evaluate solution in points x
    plot(x,y1(1,:),'b');
    hold on;
    plot(x,y1(2,:),'r');

    xlabel('*time*');
    ylabel('y in blue and dydt in red');
return;

%****** ODE y''= - (k/m)*y  ********
% nondamped oscillation
% define boundary conditions y(0)=1 and y(1)=1.5 
% writen in the form g(y(a),y(b))=0, residual
function res = bvp_bc_def(ya,yb)
    res = [ ya(1) - 1	
	    yb(1) - 1.5 ];
return;

% define ODE as a 2-D 1-st order system:
function dydt = bvp_def(t,y)
    m= 1; k=100;
    dydt = [ y(2)	
        - (k/m)*y(1)];
return;