function pensimulate(m,c,g,l,T,x0,v0,f)
%Try for example:
%pensimulate(1,0.22,1,1,1e2,0,0,@(t) 1.150*cos(2*t./3))
%pensimulate(1,0.22,1,1,1e3,0,0,@(t) 1.150*cos(2*t./3))
%pensimulate(1,0.22,1,1,1e4,0,0,@(t) 1.150*cos(2*t./3))
%(chaotic system emulates Browian motion)

G = g/l;
C = c/m;
F = @(t) f(t)/m/l;

[te,Xe] = ode15s(@ (t,y) penodefun(t,y,G,C,F),[0 T],[x0 v0]);
Xh=Xe(:,1);

figure
plot(te,Xh,'-bo');
xlabel('Time, t');
ylabel('x(t)');
title(['Numerical solution: m = ' num2str(m) ', c = ' ...
       num2str(c) ', l = ' num2str(l)]);


function z = penodefun(t,y,G,C,F)
% m*L*theta''+c*L*theta'+m*g*sin(theta)=f
% implies
% theta'=v
% v'=-c/m v -g/L*sin(theta)+f/mL
% where C=c/m, G=g/L, F=f/mL

  z=zeros(2,1);
  z(1) = y(2);
  z(2) = -C*y(2) - G*sin(y(1)) + F(t);
  
