Hello, Could anyone help me in rectifying the error while solving the ODE using bvp4c?

function main

global yright

figure (1);

hold on;

yright=0.0;

P=10e-9;

E=102e9;

h=200e-9;

b=2*h;

Es=7.56*1;

A=b*h;

v=0.34;

G=38.05e9;

I=(b*h.^3)/12;

e31=-17.05*0;

e31s=(-3e-8)*0;

k33=1.76e-8;

V=-0.2*0;

EIeff=(E+(e31.^2)/k33)*I+((1/2)*(Es+e31*e31s/k33)*b*h^2)+((1/6)*(Es+e31*e31s/k33)*h^3);

alpha=(1/((G*A*EIeff)^0.5));

J =(P^2)*alpha;

To=1*1;

H=(2*b*To)/(EIeff);

T=((e31*V*b)+(2*b*V*e31s/h))/(EIeff);

F=P/(EIeff);% here F is equal to force/(E*I)

L=9.9900E-07;

solinit=bvpinit(linspace(0,L,1000),[0 0 0 0]);

options = bvpset('NMax',1000);

sol1=bvp4c(@(x,y)cantileverode(x,y,H,T,F,L),@(ya,yb)cantileverbc(ya,yb,H,T,F,L),solinit,options);

xint=linspace(0,L,100);

Sxint1=deval(sol1,xint);

axis(['auto']);

plot(xint,Sxint1(1,:));

h = findobj(gca,'Type','line');

x=get(h,'Xdata');

y=get(h,'Ydata');

Dx=trapz(x,cos(y));

Dy=trapz(x,sin(y));

Mo=(H+T)*EIeff*(y(end)*x(end)-trapz(x,y));

U2=(J*Dx^2)+(alpha*Mo^2);

disp(Dy);

disp(U2);

xlabel('Normalized arc length (s/a)')

ylabel('Normalized rotation angle-(Phi/(Fa^2/2EI))')

function dy=cantileverode(x,y,H,T,F,L)

global yright

dy=[y(2)

   (-F*cos(y(1))-(H+T)*y(1)+(H+T)*yright)];

function res=cantileverbc(ya,yb,H,T,F,L)

global yright

yright = yb(1);

res=[ya(1)

     yb(2)];

and the error is shown to be as below;

Error using bvparguments (line 108) Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT,OPTIONS): The derivative function ODEFUN should return a column vector of length 4.

Error in bvp4c (line 129) [n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] = ...

Error in main (line 29) sol1=bvp4c(@(x,y)cantileverode(x,y,H,T,F,L),@(ya,yb)cantileverbc(ya,yb,H,T,F,L),solinit,options);

Thanks in advance

[Jan]

Please leran how you use the "{} Code" button to format your code in teh forum. Currently it is hard to read. Thanks.

Answer 1

[Jan]

Your initial solution has 4 elements:

solinit = bvpinit(linspace(0, L, 1000),[0 0 0 0]);

cantileverbc replies a [2 x 1] vector only.

[Richardson Joseph]

Thanks Jan Simon