function okrajova_uloha
t=0:0.01:1;
solinit = bvpinit(t,[0.2 0]);  % pociatocny odhad
options=[];
sol=bvp4c(@system_of_ode,@boundary_conditions,solinit,options);
solution=deval(sol,t);
plot(t,solution(1,:),'-ro',t,solution(2,:),'-.b')
legend('x','psi')

function dxpsi = system_of_ode(t,xpsi)
  x     = xpsi(1);
  psi   = xpsi(2);
  dxpsi = [ -x^3 + psi        % diferencialna rovnica pre x
             x + 3*x^2*psi ]; % diferencialna rovnica pre psi
   
function result=boundary_conditions(x0psi0,xTpsiT)
   result = [ x0psi0(1)-5     % pociat. podm. v tvare x(0)-5=0
              xTpsiT(2) ];    % koncova podm. v tvare psi(T)=0
   
