how can I solve a system of ODEs with some terminal conditions in matlab?

Question

I have 14 first order differential equations. 14 conditions, 7 are initial ones like x1(0)=0, x2(0)=5... 7 are terminal ones x8(10)=25,x9(10)=0....

I think I should use bvp4c

I found this answer but I'm still have a couple of problems:how to solve a system of Ordinary Differential Equations (ODE's) in Matlab

I create a matlab function to put my system in.

   x'=2x
   y'=3x+5y

coding it like:

 xdot=[2x(1);3x(1)+5x(2)]

like I would do in ode45. Then I should do the same for the boundary conditions. But I have no idea on how to code them. I should build a matrix containing them, but I haven't understood how to build it.

I'm trying to use this reference: http://www.math.tamu.edu/~phoward/m442/matode.pdf pag 12, but he does the y2=y' thing and I'm kind lost in my case. also it doesn't explain well how should I put the 14 conditions I have. 7 on one line and 7 on the other? how do I tell the program which value is referred to each variable?

thanks in advance.

here's the actual system. it's a bit huge, so I fear I need numerical methods.

f1=(delta1*gn-(beta*phi*x(7)*x(1)+(1-u1))/(x(1)+x(2)+x(3)+x(4))-mu*x(1)+psi*x(4));
f2=((beta*phi*x(7)*x(1)*(1-u1))/(x(1)+x(2)+x(3)+x(4))-d*x(2)-mu*x(2));
f3=(d*x(2)-(r+r0*u2)*x(3)-(alfa+mu)*x(3));
f4=((r+r0*u2)*x(3)-(mu+phi)*x(4));
f5=(delta2*hp-(phi*teta*x(3)*x(5)*(1-u1))/(x(1)+x(2)+x(3)+x(4))-gamma*x(5));
f6=((phi*teta*x(3)*x(5)*(1-u1))/(x(1)+x(2)+x(3)+x(4))-gamma*x(6)-k*x(6));
f7=(k*x(6)-gamma*x(7));
f8=x(8)*(mu - (beta*phi*x(1)*x(7) - u1 + 1)/(x(1) + x(2) + x(3) + x(4))^2 + (beta*phi*x(7))/(x(1) + x(2) + x(3) + x(4))) + x(9)*((beta*phi*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4)) - (beta*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) + (teta*x(12)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 - (teta*x(13)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f9=x(9)*(d + mu - (beta*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) - d*x(10) - A1 - (x(8)*(beta*phi*x(1)*x(7) - u1 + 1))/(x(1) + x(2) + x(3) + x(4))^2 + (teta*x(12)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 - (teta*x(13)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f10= x(10)*(alfa + mu + r + r0*u2) - A2 - x(11)*(r + r0*u2) - x(12)*((teta*phi*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4)) - (teta*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) + x(13)*((teta*phi*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4)) - (teta*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) - (x(8)*(beta*phi*x(1)*x(7) - u1 + 1))/(x(1) + x(2) + x(3) + x(4))^2 - (beta*x(9)*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f11=x(11)*(mu + phi) - x(8)*(psi + (beta*phi*x(1)*x(7) - u1 + 1)/(x(1) + x(2) + x(3) + x(4))^2) - (beta*x(9)*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 + (teta*x(12)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 - (teta*x(13)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f12=x(12)*(gamma - (teta*phi*x(3)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))) + (teta*x(13)*phi*x(3)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4));
f13=x(13)*(gamma + k) - k*x(14);
f14=gamma*x(14) + (beta*x(8)*phi*x(1))/(x(1) + x(2) + x(3) + x(4)) + (beta*x(9)*phi*x(1)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4));

extra:

u1=max(a1,min(b1,1/(2*B1)*(beta*phi/(x(1)+x(2)+x(3)+x(4))*x(7)*x(1)*(x(9)-x(8))+phi*teta/(x(1)+x(2)+x(3)+x(4))*x(3)*x(5)*(x(13)-x(12)))));
u2=max(a2,min(b2,1/(2*B2)*(r0*x(3)*x(10)-r0*x(3)*x(11))));

Answer 1

To solve BVP ode's using syms, the ode is y''+3 y' + 3 y = 0, it is first converted to 2 first order (state space formulation) and then solved

clear all; close all
syms x(t) y(t)
Dx  = diff(x);
Dy  = diff(y);
eq1 = Dx == y;
eq2 = Dy == -3*x-5*y;
[x,y] = dsolve(eq1,eq2, x(0) == 0, y(1) ==1)
figure;
ezplot(x,[0,6])

To solve same BVP using bvp4c

clear all
t0 = 0; %initial time
tf = 6; %final time
odefun=@(t,y) [y(2); -3*y(1)-5*y(2)];
bcfun=@(yleft,yright) [yleft(1);yright(1)-1];  
solinit = bvpinit(linspace(t0,tf),[0 1]);

sol = bvp4c(odefun,bcfun,solinit);

figure;
plot(sol.x(1,:),-sol.y(1,:),'r')
title('solution');
xlabel('time');
ylabel('y(t)');
grid;

ps. the numerical solution y-axis value ticks seems not to match the syms. But It looks this is just scaling of value thing. No time to look into it. May be someone can spot something and I'll update.