That is quite a long list.  You will be an expert user if you complete all the tasks and understand what you did.  Good luck.

1. Compute the paraxial ABCD matrix between 2 surfaces.

See this discussion with sample code:  https://community.zemax.com/code-exchange-10/api-python-and-matlab-interactive-extension-calculate-the-ray-transfer-abcd-matrix-for-thick-lens-1969
You might need to become an expert on this topic to truly stand by your work.  Read a bunch of textbooks on computing ABCD.
 

2. Bend a lens. Given one radius, compute the other radius to maintain the focal length.

This should be a straightforward easy first-order calculation, and you should program it in ZPL.  It is good practice for you.
 

3. Macro to plot and optionally list the distortion over the full image format

ZOS has this functionality built-in using the Grid Distortion analysis feature.

4. Calculate relative or absolute dn/dT values based on the Schott equations

ZOS has this functionality built-in using INDX MFE operand and TEMP MCO operand.  Alternatively, this should be a straightforward easy calculation, and you should program it in ZPL.  It is good practice for you.

5. Draw the aperture(s) associated with a surface.

ZOS has this functionality built-in using the Footprint analysis feature.

6. Define a macro function to compute the clearance between two surfaces, including the sags of the two surfaces

See #8 below.

7. Compute the angular boresight error in object space of the system due to the decentration of a surface or a group of surfaces.

Well, this is a bit nebulously defined, because what is your reference axis?  Boresight error in ZOS is defined in image space, with the object axis being reference.  Since you want it in object space, you are really asking to flip the system and trace from image to object.  Then you can use the BSER MFE operand.  If you do not want to flip the system, then you will need to do a bit of work and raytracing and possibly optimizing to find the chief ray angle in object space that makes it come out parallel to the axis in image space.

8. Compute the edge thickness between two surfaces at a specified x,y location

As already noted in a prior reply, Sandrine posted two separate solutions.

9. Plot the FOV points

ZOS has this functionality built-in using the Field Data Editor feature.

10. Tabulate the focal lengths of each component in the specified surface range. A component is a single surface mirror, a lens, or a cemented lens

This should be a straightforward easy first-order calculation, and you should program it in ZPL.  It is good practice for you.
 

11. Compute the local focal length in the x and y directions over the FOV.

This should be a straightforward easy first-order calculation, and you should program it in ZPL.  It is good practice for you, especially if you have done Task #10 above.
 

12. List and plot focal length vs. wavelength or the percent change in focal length vs. wavelength

ZOS has nearly this functionality built-in using the Chromatic Focal Shift analysis feature.  If you absolutely need to have the actual EFL, you will need to program it in ZPL.  It is good practice for you.

13. Do a ghost analysis between two specified surfaces

ZOS has this functionality built-in using the Ghost Focus Generator analysis feature.

14. Define a macro function to compute the size of a one-bounce or two-bounce ghost image at a specified image surface

ZOS has nearly this functionality built-in using the Ghost Focus Generator analysis feature.  If you need to customize the results, you will need to program it in ZPL.  You will need to do both forward and backward ray tracing using the paraxial raytracing equations.  It is good practice for you.  Oh, you want to optimize on it as a macro operand in the MFE.  Yes, you will need to make your own.  Here are two examples of a single bounce and double bounce ZPL.  To debug either, put in print statements all over and run it interactively, rather than as an MFE macro operand:

!zpl18.zpl

!single bounce ghost constraint macro for optimizing a problem surface

!

!g nadorff  28-jun-2015

paraxial on

! passed values:

! pvhx() == beam diameter entering system

! pvhy() == start forward raytrace at this surface

! pvpx() == reflection surface

! pvpy() == end surface to eval

!trace marginal ray:

!   RAYTRACEX x, y,        z, l, m, n, surf,   wavelength

    raytracex 0, pvhx()/2, 0, 0, 0, 1, pvhy(), 1

    s = pvpx()

    eval = pvpy()

    y = rayy(s)

!   reflect off s:

    u = -raym(s-1)/rayn(s-1)-2*y*curv(s)

!   trace backwards:

    for r = s-1, eval, -1

        y = y-thic(r)*u

        u = indx(r)/indx(r-1)*u + y*(indx(r)-indx(r-1))/indx(r-1)*curv(r)

    next

!   spot diam:

    optreturn 0 = 2*y

!   pim location after eval surface:

    optreturn 1 = y/u

!print 2*y, y/u

paraxial off

!zpl17.zpl

!double bounce ghost constraint macro for optimizing

!constraint is to achieve RMS dia/airy dia ratio > 50 after paraxial lens

!call it f = 100 mm

!

!

!Hx = first reflection surface i

!Hy = second reflection surface j

!Px = final desired surface k

!Py = unused

!

!g nadorff  7-jan-03

paraxial on

!trace marginal ray for entire system:

 raytrace 0, 0, 0, 1, 1

i = pvhx()

j = pvhy()

k = pvpx()

!   reflect off i:

    y = rayy(i)

    u = -raym(i-1)/rayn(i-1)-2*y*curv(i)

!   trace backwards from surf i to surf j+1:

    if ((i-j)==1)

       goto 1

    endif

    for s = i-1, j+1, -1

        y = y-thic(s)*u

        u = indx(s)/indx(s-1)*u + y*(indx(s)-indx(s-1))/indx(s-1)*curv(s)

    next

label 1

!   reflect off surface j

    y = y - thic(j)*u

    u = -u-2*y*curv(j)

!   trace forward to rest of surfaces

    for s = j+1, k, 1

        y = y+thic(s-1)*u

        u = indx(s-1)/indx(s)*u + y*(indx(s-1)-indx(s))/indx(s)*curv(s)

!        focus = -y/u

!        optreturn s = abso(focus)

!        optreturn s+20 = abso(y)

    next

!

!   trace thru paraxial perfect lens f=100 mm:

    u = u - y/100

    y = y+100*u

    optreturn 1 = abso(y)

paraxial off

wonderful!