So that you might comprehend how I derived the formula for d2jaMandelbrot I'm including the Maxima Cas code that I used . file iterdiff2a.wxm . /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ /* [ Created with wxMaxima version 11.08.0 ] */ /* [wxMaxima: comment start ] . (c) Copyright 2017 , sciwise@ihug.co.nz Methodology for expressing differential equations in terms of an iterated function. In particular a first order differential equation and the derivative thereof . Note that f(z,c) might be any iterated function . . [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] ' Differential equation we're examining . diff(f(z),z) + 3*f(z)*z = sin(z) for an iterated function z = f(z) , g(z) =f(f(z)) diff(g(z),z) + 3*g(z)*f(z) = sin(f(z)) . [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ kill(all)$ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ expand(eq1); expand(diff(eq1,z)) ; /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] For the iterated Mandelbrot fractal function . [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ f(z,c):=z*z+c$ g(z,c):=f(f(z,c),c)$ f(z,c); g(z,c); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ ed:expand(eq1)$ edp:expand(diff(eq1,z))$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ ed; edp; /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] This may well simplify the construction of fractint formulas for this type of equation . [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] z=Pixel c=Pixel : z=z*z+c ed = -sin(z^2+c)+3*z^6+9*c*z^4+4*z^3+9*c^2*z^2+3*c*z^2+4*c*z+3*c^3+3*c^2 edp = -2*z*cos(z^2+c)+18*z^5+36*c*z^3+12*z^2+18*c^2*z+6*c*z+4*c [wxMaxima: comment end ] */ /* Maxima can't load/batch files which end with a comment! */ "Created with wxMaxima"$ And file iterdiff2a.mac . /* . (c) Copyright 2017 , sciwise@ihug.co.nz Methodology for expressing differential equations in terms of an iterated function. In particular a first order differential equation and the derivative thereof . Note that f(z,c) might be any iterated function . . */ /* ' Differential equation we're examining . diff(f(z),z) + 3*f(z)*z = sin(z) for an iterated function z = f(z) , g(z) =f(f(z)) diff(g(z),z) + 3*g(z)*f(z) = sin(f(z)) . */ kill(all)$ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ expand(eq1); expand(diff(eq1,z)) ; /* For the iterated Mandelbrot fractal function . */ f(z,c):=z*z+c$ g(z,c):=f(f(z,c),c)$ f(z,c); g(z,c); /* */ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ ed:expand(eq1)$ edp:expand(diff(eq1,z))$ /* */ ed; edp; /* This may well simplify the construction of fractint formulas for this type of equation . */ /* z=Pixel c=Pixel : z=z*z+c ed = -sin(z^2+c)+3*z^6+9*c*z^4+4*z^3+9*c^2*z^2+3*c*z^2+4*c*z+3*c^3+3*c^2 edp = -2*z*cos(z^2+c)+18*z^5+36*c*z^3+12*z^2+18*c^2*z+6*c*z+4*c */