function comp_pendfric_lin
% Comparing solving pendfric_lin using euler with h=0.4 and ode23s

h=0.4;

figure;
hold on;    
xlim([-10 10]);
ylim([-3 3]);
xlabel('x_1');
xlabel('x_2');
title('Euler');
plot(pi*(-3:3), zeros(1,7), 'b*', 'DisplayName', 'Steady States');

x0=[-4; 1];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=eul(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

x0=[2; 2];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=eul(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

x0=[2; 1];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=eul(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

x0=[-4; 0];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=eul(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

% Comparing solving pendfric_lin using rk_classical with h=0.4 and ode23s
figure;
hold on;    
xlim([-10 10]);
ylim([-3 3]);
xlabel('x_1');
xlabel('x_2');
title('Runge-Kutta');
plot(pi*(-3:3), zeros(1,7), 'b*', 'DisplayName', 'Steady States');

x0=[-4; 1];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=rk_classical(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

x0=[2; 2];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=rk_classical(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

x0=[2; 1];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=rk_classical(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');

x0=[-4; 0];
[~,x]=ode23s(@pendfric_lin, [0 50], x0);
plot(x(:,1), x(:,2), '-r*');
[~,x]=rk_classical(@pendfric_lin, 0, 70, x0, h);
plot(x(:,1), x(:,2), '-k.');
end