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.

[time,x,y]=ex6(k,v0x,v0y)
```function [time,x,y]=ex6(k,v0x,v0y)
% Shampine's Example ch2ex3.m from
% Solving ODEs with Matlab, Cambridge University Press, 2003
% L.F. Shampine, I. Gladwell, S. Thompson, pp. 98-102, Example 2.3.5
% use [time,x,y]=ex6(0.32,0,0);
%
fill([0 0 1],[0 1 0],[0.8 0.8 0.8]);
axis([-0.1 1.1 0 4]);
hold on
plot(0,2,'ro');
n=64;dom=[0 1];kind=2;g=9.81;tol=1.e-13;
b=zeros(n,2);b(1,2)=-g;
x0=0;xp0=v0x;y0=2;yp0=v0y;%viteza initiala
time=[];x=[];y=[];
%tic;
for j=1:50
D=deriv(n,dom);T=cpv(n,dom(1),dom);t=pd(n,dom,kind);
A=D^2;A(n-1,:)=T;A(n,:)=T*D;b(n-1:n,:)=[x0,y0;xp0,yp0];
sol=A\b;solnum=t2x(sol,kind);dsol=D*sol;dsolnum=t2x(dsol,kind);
J=find(solnum(:,1)<=1);
% find the touch time tk
c=sol(:,1)+sol(:,2)-[2;zeros(n-1,1)];
[z,fz,m]=chebroots(c,dom,dom,tol);
z=real(z);z=z(abs(fz)<tol);I=length(z);
if I==0
plot(solnum(J,1),solnum(J,2),'r');hold on;%pause(0.2);
x0=solnum(end,1);xp0=dsolnum(end,1);
y0=solnum(end,2);yp0=dsolnum(end,2);
else tk=z(1);
Tr=cpv(n,tk,dom);
x0=Tr*sol(:,1);xp0=-k*Tr*(D*sol(:,2));
y0=Tr*sol(:,2);yp0=k*Tr*(D*sol(:,1));
J1=find(t<tk);J2=intersect(J,J1);
if x0<=1,
plot([solnum(J2,1);x0],[solnum(J2,2);y0],'r',x0,y0,'og');
hold on;%pause(0.2);
else
plot(solnum(J2,1),solnum(J2,2),'r');hold on;%pause(0.2);
end
if x0>=1, break;end
if ~isempty(time)&&(tk-time(end))/time(end)<0.001,
display('Cluster, end bouncing');break;
end
time=[time;tk];x=[x;x0];y=[y;y0];dom=[tk,tk+1];
end
end
%toc;
grid;xlabel('x');ylabel('y');hold off;
display([time,x,y]);```