Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

15 Jul 2011 (Updated )

The MATLAB package Chebpack solves specific problems for differential or integral equations.

[xg,solnumg]=ibvp_sys_succapprox_ex1(n,dom,kind,y0)
```function [xg,solnumg]=ibvp_sys_succapprox_ex1(n,dom,kind,y0)
% Example: y1'=B-(A+1)y1+y1^2 y2, y2'=Ay1-y1^2 y2, x in [0,200]
%          y1(0)=1.5, y2(0)=3, A=1.9, B=1
%          (the Brusselator)
% test: fixed point at (B,A/B)
% call: [x,solnum]=ibvp_sys_succapprox_ex1(64,[0,1,2:2:200],2,[1.5,3]);
%
tic;xg=[];solnumg=[];solg=[];p=length(dom)-1;dyg=[];contg=0;m=length(y0);
A=spalloc(m*n,m*n,10*p*m);
for i=1:p
domloc=[dom(i),dom(i+1)];
x=pd(n,domloc,kind);X=mult(n,domloc);J=prim(n,domloc);xg=[xg;x];
myDElin;y=zeros(m*n,1);cont=1;dy=1;
for k=1:m,y((k-1)*n+1)=2*y0(k);end
while dy(cont) > 1.e-12
if cont>100, break, end;
cont=cont+1;myDEnonlin;sol=A\b;dy(cont)=norm(sol-y);y=sol;
end
solnum = t2x(reshape(sol,n,m),kind);
solnumg=[solnumg;solnum];solg=[solg;sol];dyg=[dyg,dy];contg=contg+cont;
y0=solnum(end,:);
end;toc;myOUT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDElin
% describes the linear part of the differential system
% implements the initial conditions in the matrix A
A=[speye(n)+2.9*J,zeros(n);-1.9*J,speye(n)];
T=cpv(n,domloc(1),domloc);A(1,:)=[T,zeros(1,n)];A(n+1,:)=[zeros(1,n),T];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDEnonlin
% describes the nonlinear part of the differential system
% implements the initial conditions in the vector b
yv=t2x(reshape(y,n,2),kind);work=(yv(:,1).^2).*yv(:,2);
b=[J,zeros(n);zeros(n),J]*reshape(x2t([1+work,-work],kind),2*n,1);
b(1)=y0(1);b(n+1)=y0(2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
% describes the output of the code
figure(1);
subplot(2,1,1);semilogy(abs(solg),'.');
title('Absolute value of the coefficients of y_1, y_2');grid;
subplot(2,1,2);plot(xg,solnumg(:,1),'b',xg,solnumg(:,2),'r--');
title('Numerical solution');legend('y1(x)','y2(x)');
xlabel('x');ylabel('y1,2(x)');grid;
figure(2);
semilogy(1:contg,dyg,'.');grid;title('History of iterations - corrections');
figure(3);
plot(solnumg(:,1),solnumg(:,2),'b',1,1.9,'.r');grid;
xlabel('y1(x)');ylabel('y2(x)');title('Phase portrait');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end```