function sol=bvp4cfd(odefun,bcfun,solinit,options) if (isfield(solinit,"x")) if columns(solinit.x)<5 error("bvp4c: initial mesh solinit.x must have 5 or more points."); endif t = solinit.x; else error("bvp4c: missing initial mesh solinit.x"); end if (isfield(solinit,"y")) y = solinit.y; else error("bvp4c: missing initial guess"); end if (isfield(solinit,"parameters")) error("bvp4c: solving for unknown parameters is not yet supported"); end RelTol = 1e-3; AbsTol = 1e-6; if ( nargin > 3 ) if (isfield(options,"RelTol")) RelTol = options.RelTol; endif if (isfield(options,"RelTol")) AbsTol = options.AbsTol; endif endif ##profile off; ##profile clear; ##tic; ##profile on; nvar=size(y,1); maxref=14; for kk=1:maxref h=diff(t); np=length(t); ci=(kron(1:5,ones(np,1))+min(floor([0 0 0:np-3]'),np-5)).'; c=zeros(5,np); ch=zeros(5,np); dch=zeros(5,np); tcit=(t(ci)-(t)).'; [utcit,i,j]=unique(tcit,"rows"); for k=1:rows(utcit) kj=k==j; sumkj=sum(kj); c(:,kj)=kron(finitedifferencecoef1d(utcit(k,:),1),ones(1,sumkj)); if i(k)=AbsTol)&&(RelErr>=RelTol))&&kk=AbsTol)&(err./max(max(abs(dy)))>=RelTol); t_add=tm(ref_int); y_add=ym(:,ref_int); [t ti]=sort([t,t_add]); y=[y y_add](:,ti); else break end end sol.solver="bvp4c"; sol.x=t; sol.y=y; sol.yp=dy; ##profile off; ##toc ##profshow () ##profexplore () end function r=res(y,ci,C,nvar,t,odefun,bcfun) y=reshape(y,nvar,length(t)); r=[(y(:,ci)*C-odefun(t,y))(:);bcfun(y(:,1),y(:,end))]; end function C=finitedifferencecoef1d(c,m) %Calculates the coefficients C_i of the weighted sum to obtain the finite difference approximation of the derivative % sum_(i=0)^(n-1)C_i f(x-c_i) %where c_i is the distances from x %Example input to C=finitedifferencecoef1d(c,m) %m=8; %c=linspace(-2,2,20); %Inputs %c is the 1 by n vector of distances from x %m is the order of the derivative % %An example of using the result to obtain the mth derivative: %x=.1; %sin(x+c)*C %A comparison with the true derivative %t=[sin(x) % cos(x) % -sin(x) % -cos(x)]; %t(rem(m,4)+1) n=size(c,2); M=zeros(n,1); M(m+1)=round(gamma(m+1));%1;%factorial(m) S=(kron(c,ones(n,1))).^kron((0:n-1)',ones(1,n));%./repmat(factorial((0:n-1)'),1,n); C=S\M; end %!demo %! t=linspace(0,4,5); %! ic=bvpinit(t,[0;0]); %! dydt=@(t,y)[y(2,:);-abs(y(1,:))]; %! bcfun=@(ya,yb)[ya(1);yb(1)+2]; %! s=bvp4cfd(dydt,bcfun,ic); %! [tt,yy]=ode45(dydt,[0 4],s.y(:,1)); %! plot(s.x,s.y,'-d',tt,yy,'-x'); %!demo %! t=linspace(1/3/pi,1,5); %! ic=bvpinit(t,[0;0]); %! dydt=@(t,y)[y(2,:);-y(1,:)./t.^4-2./t.*y(2,:)]; %! bcfun=@(ya,yb)[ya(1);yb(1)-sin(1)];; %! s=bvp4cfd(dydt,bcfun,ic); %! tt=s.x; %! yy=[sin(1./tt);-cos(1./tt)./tt.^2]; %! disp("maximum error"); %! disp(max(abs(s.y-yy),[],2)); %! [ttt,yyy]=ode45(dydt,[1 1/(3*pi)],s.y(:,end)); %! plot(s.x,s.y,'-d',tt,yy,'-',ttt,yyy,'-o');