First derivative can be calculated by the complex-step derivative formula:
$f'(x)=\frac{Im(f(x+ih))}{h}$
Generalization of the above for calculating derivatives of any order employs multicomplex numbers, resulting in multicomplex derivatives:
$f^{(n)}(x)=\frac{C^{(n)}_{n^2-1}(f(x+i^{(1)}h+i^{(n)}h))}{h^n}$
According to the Wiki Complex-variable methods:
In Matlab, the calculation of the first order derivative is very easy to implement:
x=0:0.01:10;
h=0.001;
f=sin(x);
df=imag(sin(x+h*i))./h;
plot(x,f)
hold on
plot(x,df)
I do not understand how to implement the calculation of the second order derivative, because I do not understand what is $i^{(1)},i^{(2)}...i^{(n)}$ and how the operator $C^{(n)}_{n^2-1}$ is calculated.
EDIT: Here is the program for the second order derivative, which is calculated incorrectly. I don't understand how to use $Imag_{12}$
x=0:0.01:15;
h=0.0000001;
imx1=x+(i)*h;
imx2=x+(i+i)*h;
f=(x).^2;
df = imag((imx1).^2)./h;
ddf = imag((imx2).^2)./h^2;
$ \def\bbR#1{{\mathbb R}^{#1}} \def\LR#1{\left(#1\right)} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} $Since you're using Matlab, a simpler approach is to create a Jordan block matrix $$\eqalign{ &J = (xI + N) \; \in\,\bbR{n\times n} \\ &N^n = 0 \qquad \big\{{\rm nilpotent}\big\} }$$ Have Matlab evaluate the function using this matrix as its argument, and since $$f(J) = \sum_{k=0}^{n-1}\frac{N^k\,f^{(k)}(x)}{k!}$$ the $k^{th}$ derivative can be read off of the $k^{th}$ superdiagonal of the result for $\,1\le k<n$