function [j]=jacob(x,u,deltax,parameter)
%function [j]=jacob(x,u,deltax)
% j is the jacobian of the functions in function fun.m
% x in the state vector in fun(x,u,parameter)
% u is a vector of constants
% j is square the size of x
[n,m]=size(x);
j=zeros(n);
y0=fun(x,u,parameter);
for i=1:n
	x1=x;
	x1(i)=x1(i) + deltax(i);
	y1=fun(x1,u,parameter);
	dy=y1 - y0;
	dy=dy/deltax(i);
	j(:,i)=dy;
end
