working with stellar spectra in fits format in python
The information you need to recreate the wavelength array is in the World Coordinate System (WCS) of the header, specifically:
CRPIX1  =                 1.00
CRVAL1  =            3500.0000 / central wavelength of first pixel
CDELT1  =             0.900000 / linear dispersion (Angstrom/pixel)

which lists the starting/reference pixel of the wavelength array (1.0), the wavelength value at the start point (3500 angstroms (assumed)) and the step per pixel (0.9 Angstrom/pixel). To read this information, it is best to use a WCS library rather than trying to interpret them directly as they can be more complicated and there are many subformats of FITS WCS.
Fortunately astropy has a module to make this easy (starting from your code above and extending it):
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS

hdul = fits.open('s0013.fits')
data = hdul[0].data
h1 = hdul[0].header
obj_name = h1.get('OBJECT', 'Unknown')

flux = data[0]
w = WCS(h1, naxis=1, relax=False, fix=False)
lam = w.wcs_pix2world(np.arange(len(flux)), 0)[0]

plt.plot(lam, flux)
plt.ylim(0, )
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Normalized flux')
plt.savefig(obj_name + '.png')

This will produce the following plot: 
If you want to do more extended manipulation of spectra, particularly with the bewildering variety of wavelength and flux units, it might be worth looking at synphot and specutils which build on Astropy and add more direct support for spectra beyond simple numpy arrays. For example, you could make a synphot SourceSpectrum from the above by doing:
from astropy import units as u
from synphot import units, SourceSpectrum
from synphot.spectrum import Empirical1D

source_spec = SourceSpectrum(Empirical1D, points=lam*u.AA, lookup_table=flux,
                   keep_neg=True, meta={'header': h1})