Choose better initial conditions e.g. for q1 and q2 they are both zero
close all
D=0.02;
r=0.04;
A=2.5;
B=0.01;
z=3.0;
C=3.25;
%>>>>>>>>> choose hetter intial coditions
q1=1;
q2=0.5;
%-->>>>>>>>>
p1=0.632;
E=0.2;
p2=sqrt(2*E- p1^2);
dt=1/5;
tspan=[10:dt:100];
K=[0 0.050 0.5 1.5];
m=length(K);
figure(1)
hold on
for i=1:m;
[t,x]=ode45(@(t,x)IP(t,x,B,C,D,r,z,A,K(i)),[0 100],[q1 q2 p1 p2]);
subplot(2,2,(i))
plot(x(:,1),x(:,3),'-b',x(:,2),x(:,4),'-r ')
end
xlabel('q');
ylabel('p');
function dxdt = IP(t,x,B,C,D,r,z,A,K)
dxdt = zeros(4,1);
dxdt(1) = x(3);
dxdt(2) = x(4);
dxdt(3) = (z*D*exp(-z*x(1))*(2*A*exp(2*z*r)*exp(-z*x(1)))+(B*exp(z*r))*(1+exp(-(z*x(1))))+2*C)/((1-exp(-(z*x(1))))^3)+K*(x(2)-x(1));
dxdt(4) = (z*D*exp(-z*x(2))*(2*A*exp(2*z*r)*exp(-z*x(2)))+(B*exp(z*r))*(1+exp(-(z*x(2))))+2*C)/((1-exp(-(z*x(2))))^3)+K*(x(1)-x(2));
end
