You have the gist of it, but there are a couple of problems. I'm sorry, I don't have MATLAB, but I'll do my best to give you MATLAB-like pseudo-code to work with.
The main problem is that your Airy disk is an image, i.e. amplitude only (the result of taking `abs`). To reconstruct the circular aperture, you need phase as well. Here's code to first generate your complex Airy disk (PSF), from which you can then reconstruct the circular aperture (pupil function). Starting with this file, `circ.jpg`:
pupil = imread('circ.jpg');
pupil = rgb2gray(pupil);
% use fftshift to block-swap the fft output (move DC to the center)
complex_airy = fftshift(fft2(pupil));
% now, starting with this complex airy pattern, you were close:
pupil_reconstructed = ifft2(complex_airy);
pupil_reconstructed = abs(pupil_reconstructed);
Below are the amplitude and phase of the complex_airy pattern, obtained using `abs(complex_airy)` and `angle(complex_airy)`, respectively. The problem with your initial approach is that it's an ill-formed inverse problem: a single Airy amplitude pattern can be generated from an infinite number of pupil functions. Of course there are statistical approaches for estimating the pupil function from the Airy amplitude, using lots of other constraints gleaned from the light source properties and optical system, but I don't think that's what you were after.
