0001 function [ final ] = stencil_order( A,order )
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 [m,n,l]=size(A);
0012 if m~=n;
0013 error('input array must be square or cube');
0014 elseif mod(n,2)~=1;
0015 error('matrix must have odd side length');
0016 elseif mod(order,1)~=0;
0017 error('input order must be an integer value');
0018 end
0019
0020 c=(n+1)/2;
0021 syms Dx Dy Dz h;
0022 e=0;
0023
0024 order = order+1;
0025 if l==1;
0026 for i=1:n;
0027 for j=1:n;
0028 tay=taylor(exp((i-c)*Dy*h+(j-c)*Dx*h),Dy,'order',order);
0029 tay=taylor(tay,Dx,'order',order);
0030 e=e+A(i,j)*tay;
0031 end
0032 end
0033 elseif l==n;
0034 for i=1:n;
0035 for j=1:n;
0036 for k=1:n;
0037 tay=taylor(exp((i-c)*Dy*h+(j-c)*Dx*h+(k-c)*Dz*h),Dy,'order',order);
0038 tay=taylor(tay,Dx,'order',order);
0039 tay=taylor(tay,Dz,'order',order);
0040 e=e+A(i,j,k)*tay;
0041 end
0042 end
0043 end
0044 else
0045 error('input array must be square or cube');
0046 end
0047
0048 [cx,tx]=coeffs(e,h);
0049 cx=flip(cx);
0050 tx=flip(tx);
0051
0052 final=blanks(1);
0053 hind=1;
0054 first=1;
0055 for k=1:order
0056
0057 syms a;
0058 a=k;
0059 if hind<=length(cx);
0060 if isequaln(tx(hind),h^a);
0061 temp1=char(h^a);
0062 temp2=char(cx(hind));
0063 if first==1;
0064 final=strcat(final,temp1,'{',temp2,'}');
0065 first=0;
0066 else
0067 final=strcat(final,' +',temp1,'{',temp2,'}');
0068 end
0069 hind=hind+1;
0070 end
0071 end
0072
0073 end
0074 end
0075
0076
0077
0078