function fermi_pasta()
	
	c = 50;
	q = [	1	0	0	0	-5	0	0	0	1	]'; % 10 is fixed right border!
	p = [	0	0	0	0	0	0	0	0	0	]';
	
	function energy = hamiltonian(p, q)
		energy = 0;
		for i = 1:length(p)/2
			energy = energy + 1/2 * (p(2*i - 1)^2 + p(2*i)^2) + c*(q(2*i)-q(2*i - 1))^2;
			if 2*i < length(p)
				energy = energy + (q(2*i + 1) - q(2*i))^4;
			end
		end
		energy = energy + q(1)^4;
		energy = energy + (q(end))^4;
	end
	
	function j = jacobian(p, q)
		j = zeros(length(p) * 2, 1);
		for i = 1:length(p)
			% Derivative for t of q_i → p_i
			j(length(p) + i) = p(i);
			
			% Derivative for t of p_i → -q_i
			for h = 1:length(p)/2
				if i == 2*h
					j(i) = j(i) - c * 2 * (q(2*h) - q(2*h-1));
					if i < length(p)
						j(i) = j(i) + 4 * (q(2*h+1) - q(2*h))^3;
					end
				elseif i == 2*h - 1
					j(i) = j(i) + c * 2 * (q(2*h) - q(2*h-1));
				elseif i == 2*h + 1
					j(i) = j(i) - 4 * (q(2*h+1) - q(2*h))^3;
				end
			end
		end
		j(1) = j(1) - 4*q(1)^3;
		j(length(p)) = j(length(p)) + 4 * (- q(end))^3;
	end
	
	figure('name', 'fermi_pasta');
	subplot(2, 1, 2);
	title('Energy');
	xlabel('Time t');
	ylabel('Energy');
	xlim([0 100]);
	ylim([0 5000]);
	function status = plot_springs(ts, ys, action)
		if ~isempty(action)
			return
		end
		
		for i = 1:length(ts)
			y = ys(:, i);
			t = ts(:, i);
		
			p = y(1:length(p));
			q = y(length(p)+1:end);
			
			subplot(2, 1, 1);
			plot(interp(linspace(0, 10, length(p) + 2), 3), interp([ 0; q; 0 ], 3), '-', 'LineWidth', 2);
			title('Positions');
			hold on
			grid
			scatter(linspace(0, 10, length(p) + 2), [ 0; q; 0 ], 'MarkerFace', 'red');
			ylim([ -20 20 ]);
			xlim([ 0 10 ]);
			hold off
			drawnow;
						
			subplot(2, 1, 2);
			hold on;
			plot([ t ], [ hamiltonian(p, q) ], '.', 'MarkerSize', 1);
			hold off;
			drawnow;
		end
		
		if isempty(findobj('type', 'figure', 'name', 'fermi_pasta'))
			status = 1;
			close(findobj('type', 'figure'));
			return;
		end
		
		status = 0;
	end
	
	% Some basic integrators matlab does not have per se
	function [t, x] = explicit_euler(jacobian, time, initial, params)
		step_size = time(2) - time(1);
		for t = time
			initial = initial + step_size * jacobian(t, initial);
			params.OutputFcn(t, initial, '');
		end
	end

	% explicit_euler((@(t, x) jacobian(x(1:length(p)), x(length(p)+1:end))), 0:0.001:100, [ p; q ], odeset('OutputFcn', @plot_springs));
	ode45((@(t, x) jacobian(x(1:length(p)), x(length(p)+1:end))), 0:0.01:100, [ p; q ], odeset('OutputFcn', @plot_springs));
end

