function RamseyModel
clear all
close

subplot(2,2,1);
t=0:0.01:5;
r1=Ramsey(t,0.25,0.25,0.5,2);
plot(t,r1(1,:),'-b',t,r1(2,:),'-.b');

subplot(2,2,3);
t=0:0.01:20;
r2=Ramsey(t,0.25,0.25,0.5,2);
plot(t,r2(1,:),'-r',t,r2(2,:),'-.r');

subplot(2,2,2);
plot(r1(1,:),r1(2,:),'-b',r2(1,:),r2(2,:),'-r')

function result=Ramsey(t,r,d,k0,kT)
solinit = bvpinit(t,[0.5 1]);  %initial guess
options=[];
sol=bvp4c(@system_of_ode,@bound_cond,solinit,options,r,d,k0,kT);
result=deval(sol,t);

function dkc=system_of_ode(t,kc,r,d,k0,kT)
   dkc = zeros(2,1);     % stlpcovy vektor
   dkc(1)=sqrt(kc(1))-d*kc(1)-kc(2);
   dkc(2)=-kc(2)*(d+r-1/(2*sqrt(kc(1))));
   
function result=bound_cond(kc_0,kc_T,r,d,k0,kT)
   result=zeros(2,1);
   result(1)=kc_0(1)-k0; % pociatocna podmienka v tvare k(0)-k0=0
   result(2)=kc_T(1)-kT; % koncova podmienka v tvare k(T)-kT=0
   
