function sols = hkbsimhist (b,d)

clf;


a = 1; %Coefficients determine shape of potential...

el = -4; %determines the range of possible start values for a simulation,
%which are [el:er]
er = 4;
% f = 2;  higher frequency oscillator
T = 3; %total time
N = 1000;  %number of iterations per simulation
n = 1000; %number of times to run the simulation
h = 100; %number of histogram bins
init = 0; %initial value if not random

%plot the potential using given coeffients
xp = el:.1:er;
yp = -a*cos(xp) - b*cos(2*xp);
figure (10)
subplot(2,1,1);
plot (xp,yp);
title('potential')
ylim([el er]);
%randn('state',100) % set the state of randn - uncomment to make simulation
%same every time
dt = T/N; %time step
dX = zeros(1,N); % preallocate arrays ...
X = zeros(1,N); % for efficiency
sols = zeros(1,n);

for i = 1:n,
    X(1) = ((er-el)*rand)+el; %initial value randomly chosen
    %X(1) = init; %initial value specifically set
    for j = 2:N,
        dX(j-1) = dt*(-a*sin(X(j-1)) - b*2*sin(2*X(j-1)) + d*randn); % general increment
        X(j) = X(j-1) + dX(j-1);
    end;
    sols(i)= X(N) ;
end;

%plot histogram of simulation solutions
subplot(2,1,2);
hist(sols, h);
title('histogram of solutions')
set(gca,'XLim',[el,er]);

figure (11)
plot ([1:N], X)
title('Time fuction of Last Simulation', 'FontSize', 18)
set(gca,'FontSize', 18);
xlabel ('Iteration');
ylabel ('X');

shg;