
global C K R

C = 1;
K = 1;
R = 1/2;

atol = 1e-6;
rtol = 1e-3;
opt = odeset('RelTol', rtol, 'AbsTol', atol);
opt = odeset(opt, 'Events', @odef3_events);
opt = odeset(opt, 'Mass', @odef3_mass);
opt = odeset(opt, 'MStateDependence', 'none');
opt = odeset(opt, 'MassSingular', 'yes');

w0 = 2;
Y0 = [w0; w0/(-2*R)];

% remember this code only works if you start at 1/R
%[ETA, Y] = ode45(@odef3_rhs, [1/R  1/R+20], Y0, opt);
[ETA, Y, ETAe, Ye, Ie] = ode15s(@odef3_rhs, [1/R  -1/R+.5], Y0, opt);

zcrit1 = 0;
zcrit2 = 1/R;
zcrit3 = -1/R;


figure(1);  clf;   hold on;
xlabel('\eta');
plot(ETA, Y(:,1), 'b.-');
plot(ETA, Y(:,2), 'y.-');
%plot(ETA2, Y2(:,1), 'g.-');
%plot(ETA2, Y2(:,2), 'm.-');
ax = axis;
legend('w(z)', 'w''(z)', 0);
plot([zcrit1, zcrit1], [ax(3) ax(4)], 'ro:');
plot([zcrit2, zcrit2], [ax(3) ax(4)], 'r*:');
plot([zcrit3, zcrit3], [ax(3) ax(4)], 'r*:');
grid on;
axis([-1/R 1/R -10 10]);
if ~isempty(ETAe)
  plot(ETAe, Ye(:,1), 'r*');
  plot(ETAe, Ye(:,1), 'ro');
end