This is code that I've written in Maxima CAS . This finds the first and second derivative of an iterated function , g(z,c) . To construct a fractal substitute z:g(z,c) with the equation for g(z,c) , rather than calling a function. Similarly for z1:g1(z,z1,c) and z2:g2(z,z1,z2,c) . I've provided an example of what this should look like. The validity of this approach is tested by symbolically calculating the second derivative at each iteration and comparing this with the second derivative from the above method. This is done for three iterations and , if correct , will display TRUE. The file difffrac1h.mac /* . (c) 2016 , sciwise@ihug.co.nz . A general method for defining the first and second derivatives of an iterated fractal function , g(z,c) . Some care is required to ensure that g(0,c) is defined . . */ /* Example Fractint code . d2zmand(XAXIS){; Edward Montague ; Mandelbrot series = z ; 1st Derivative of Mandelbrot series = z1. ; 2nd Derivative of Mandelbrot series = z2. z = Pixel z1 = 1 z2 = 0: z2 = 2*z1*z1 +2*z2*z z1 = 2*z*z1 + 1 z = z*z+Pixel |z2| z1 */ eq:diff(g(z,c),z,1)*z1 + diff(g(0,c),c); define(g1(z,z1,c),eq)$ /* Second Derivative --> z2 */ eq:diff(g(z,c),z,2)*z1*z1 + diff(g(z,c),z)*z2 + diff(g(0,c),c,2); define(g2(z,z1,z2,c),eq)$ /* Initial values */ z : c; z1 : diff(z,c); z2 : diff(z,c,2); /* First iteration. */ z2 : g2(z,z1,z2,c); z1 : g1(z,z1,c)$ z : g(z,c); eq:diff(z,c,2); is(equal(eq,z2)); /* Second iteration. */ z2 : g2(z,z1,z2,c); z1 : g1(z,z1,c)$ z : g(z,c); eq:diff(z,c,2)$ is(equal(eq,z2)); /* Third iteration. */ z2 : g2(z,z1,z2,c); z1 : g1(z,z1,c)$ z : g(z,c)$ eq:diff(z,c,2)$ is(equal(eq,z2)); /* */ And the file difffrac1h.wxm /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ /* [ Created with wxMaxima version 11.08.0 ] */ /* [wxMaxima: comment start ] . (c) 2016 , sciwise@ihug.co.nz . A general method for defining the first and second derivatives of an iterated fractal function , g(z,c) . Some care is required to ensure that g(0,c) is defined . . [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] Example Fractint code . d2zmand(XAXIS){; Edward Montague ; Mandelbrot series = z ; 1st Derivative of Mandelbrot series = z1. ; 2nd Derivative of Mandelbrot series = z2. z = Pixel z1 = 1 z2 = 0: z2 = 2*z1*z1 +2*z2*z z1 = 2*z*z1 + 1 z = z*z+Pixel |z2| z1 [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ eq:diff(g(z,c),z,1)*z1 + diff(g(0,c),c); define(g1(z,z1,c),eq)$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Second Derivative --> z2 [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ eq:diff(g(z,c),z,2)*z1*z1 + diff(g(z,c),z)*z2 + diff(g(0,c),c,2); define(g2(z,z1,z2,c),eq)$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Initial values [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ z : c; z1 : diff(z,c); z2 : diff(z,c,2); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] First iteration. [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ z2 : g2(z,z1,z2,c); z1 : g1(z,z1,c)$ z : g(z,c); eq:diff(z,c,2); is(equal(eq,z2)); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Second iteration. [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ z2 : g2(z,z1,z2,c); z1 : g1(z,z1,c)$ z : g(z,c); eq:diff(z,c,2)$ is(equal(eq,z2)); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Third iteration. [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ z2 : g2(z,z1,z2,c); z1 : g1(z,z1,c)$ z : g(z,c)$ eq:diff(z,c,2)$ is(equal(eq,z2)); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* Maxima can't load/batch files which end with a comment! */ "Created with wxMaxima"$