
global C K R

C = 1;
K = .01;
R = 2;

opt = odeset('Jacobian', @odef1_jac);
opt = odeset('Mass', @odef1_mass);
% Mass matrix is only dependent on t (z) but not the state
opt = odeset(opt, 'MStateDependence', 'none');
% it is singular
opt = odeset(opt, 'MassSingular', 'yes');
opt = odeset(opt, 'RelTol', 1e-5, 'AbsTol', 1e-7);


w0 = [0; .1];

%[Z, W] = ode15s(@odef1_rhs, [51 99:.01:101 149], w0, opt);
[Z, W] = ode15s(@odef1_rhs, [51 160], w0, opt);
%[Z, W] = ode45(@odef1_rhs, [52 148], w0, opt);

zcrit1 = C/K
zcrit2 = C/K - 1/(R*K)
zcrit3 = C/K + 1/(R*K)
 

figure(1);  clf; 
subplot(2,1,1); hold on;
plot(Z, W(:,1), 'b.-');
%plot(Z, W(:,2), 'g--');
xlabel('z');
ylabel('w(z)');

ax = axis;
ymin=ax(3);  ymax=ax(4);
plot([zcrit1 zcrit1], [ymin, ymax], 'ro:');
plot([zcrit2 zcrit2], [ymin, ymax], 'r*:');
plot([zcrit3 zcrit3], [ymin, ymax], 'r*:');
ax(1) = ax(1)-(ax(2)-ax(1))/20;
ax(2) = ax(2)+(ax(2)-ax(1))/20;
axis(ax);
%legend('w(z)', 'w''(z)', 0);
%legend('w(z)', 'w''(z)', 0);

subplot(2,1,2); hold on;
%plot(Z, W(:,1), 'b.-');
plot(Z, W(:,2), 'g.--');
xlabel('z');
ylabel('w''(z)');

ax = axis;
ymin=ax(3);  ymax=ax(4);
plot([zcrit1 zcrit1], [ymin, ymax], 'ro:');
plot([zcrit2 zcrit2], [ymin, ymax], 'r*:');
plot([zcrit3 zcrit3], [ymin, ymax], 'r*:');
ax(1) = ax(1)-(ax(2)-ax(1))/20;
ax(2) = ax(2)+(ax(2)-ax(1))/20;
axis(ax);
%legend('w(z)', 'w''(z)', 0);
